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The  effort  described  in  this  report  is  a  sampling  of  the  most  interesting 
problems  undertaken  in  this  contract.  Many  of  the  problems  were  a  continua¬ 
tion  in  kind  of  work  done  under  Air  Force  Contract  F19628-73-C-0136.  That 
effort  was  summarized  in  a  Final  Report  dated  31  March  1976,  titled  Numerical 
and  Data  Analysis  Techniques  Applied  to  Scientific  Research-II,  AFGL-TR-76- 
0091.  jjltf  O  '/I 

The  authors  wish  to  express  their  thanks  to  Dr.  Paul  Tsipouras  of  the 
Air  Force  Geophysics  Laboratory  (AFGL)  for  his  assistance  in  the  development 
of  solutions  to  many  of  the  problems  considered  under  this  contract.  A 
similar  debt  of  gratitude  is  owed  to  all  Problem  Initiators  from  AFGL  for 
providing  the  detailed  background  information  regarding  their  particular 
problems . 

Mr.  Leo  F.  Power,  Jr. ,  the  Director  of  this  laboratory,  was  invaluable 
in  supplying  the  necessary  resources  to  meet  critical  deadlines  and  finally, 
our  appreciation  goes  out  to  the  support  staff  of  this  laboratory  for  their 
continual  assistance  in  the  program  preparation  of  these  problems  and  in 
particular  to  Miss  Mary  Kelly  for  the  excellent  job  of  preparing  this  report 
for  publication. 
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EARTH'S  GEOIV  AND  GRAVITY  FIELD 


IwiticutoK:  Ato..  G.  Hadgyigzoxgz 

VhobLvm  No:  4919  {^otimesiZii  VtiobZm  No:  4850) 

Rnoje-dt  No:  7 600 


Of  special  interest  to  the  Air  Force  Geophysics  Laboratory  (AFGL) 
through  the  years  has  been  the  reduction  of  satellite  measurements  to  the 
various  parameters  measured.  This  is  most  certainly  true  for  satellite 
measurements  of  the  oceanic  geoid  surface.  Of  primary  concern,  is  that 
approach  which  may  remove  the  necessity  for  extremely  accurate  reference 
orbits  in  order  to  exploit  satellite  altimetry  of  1  meter  accuracy  for 
geoid  improvement.  Conventional  approaches  utilize  long-arc  orbit  integra¬ 
tion  that  requires  extensive  tracking  from  ground  based  trackers.  Whereas, 
the  short -arc  approach ^ in  program  SARRA  has  shown  that  accuracies  of 
integration  of  better  than  1  meter  can  be  attained  provided  that 

a.  spherical  harmonics  at  least  through  (n,m)  =  (4,4)  are  exercised 
in  the  integration 

b.  all  six  orbital  parameters  (x,y,z,x,y,z)  at  mid-arc  are  free 
to  adjust. 


At  present,  there  are  two  options  for  selecting  the  mathematical 
model  defining  the  geoid  surface.  The  first  is  the  spherical  multiquadric 
that  is  based  on  least  squares  fitting  of  the  observational  data.  If  the 
area  being  studied  is  covered  with  an  adequately  dense  set  of  measurements, 
then,  this  approach  has  no  limit  as  to  the  size  of  the  area  to  be  pro¬ 
cessed  or  the  detail  to  be  obtained.  There  are  two  distinct  but  coupled 
mathematical  formulations  adhered  to  in  the  programming. 


A.  The  basic  observation  equation^-*  is 


f (H. . ,x . ,y . , z . ,x . ,y . , z . ,r  ,c  )  =  o 
ij  3  3  3  3  3  3  o  n 
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where,  H. .  is  the  altimetry  measurement 


(x. ,y . ,z . ,x. ,y . ,z.) 

3  3  3  3  3  3 

the  state  vector;  r  the  reference  radius;  cn  the  geoid  coefficient. 


B.  The  geoidal  model  There  expressed  in  terms  of  spheroidal  multi¬ 
quadric  functions) (5) 
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where,  r^  is  the  geocentric  radius;  c^j  is  to  be  determined; 
is  called  the  kernel  function  and  of  the  form 
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|^(x.-x^+ka)  +  (y.-y^+ka)  +  (z^z^+ka)  : 


xi ,yi , 2i  is  the  geoidal  point;  xj,yj,zj  is  the  node  point;  a  is  the 
semi -major  axis;  and  k  is  the  arbitrary  fraction. 


The  other  is  the  Covariance  function  whose  advantage  is  realized  when 
data  is  irregularly  distributed,  with  some  areas  being  sparse  in  measurements. 
This  function  DO)  is  obtained  by  averaging  the  product  of  undulations  sepa¬ 
rated  by  the  spherical  distance  ip.  It  can  be  expressed  as 


where. 
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corrections  to  the  spherical 


harmonic  coefficient 


R  the  mean  radius  of  the  earth 
G  the  mean  gravity  of  the  earth  at  the  surface. 
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In  an  effort  to  decrease  computational  requirements  the  programmed 
mathematical  formulation  for  the  radial  distance  r^  was  replaced  by  the 
following: 


r. 
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in  which 
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equatorial  radius  of  adopted  reference  spheroid 
gravitational  constant  of  earth 

geocentric  latitude,  longitude  and  radius  on  the  mean 
sea  level  geoid 

coefficients  of  the  potential  function  expanded  in 
spherical  harmonics 

associated  Legendre  polynomial 

rotational  velocity  of  earth. 


The  solution  of  r^  is  to  be  obtained  from  an  iterative  process  using  the 
following  for  a  starting  value: 


where  a  and  <f>  are  as  defined  previously,  and  e  represents  the  eccentricity. 

The  observed  data  was  obtained  from  in  excess  of  100  tapes  of  altimetry 
data  received  from  NASA  and  NSWC  (Naval  Surface  Weapons  Center) .  This  data 
was  reduced  and  edited  by  the  preprocessor  computer  program  which  builds 
the  input  data  base  for  the  SARRA  reductions. 

The  preprocessing  and  editing  (PREP)  of  GEOS-3  altimetry  data  is 
designed  as  a  three  level  effort.  The  raw  data  is  first  edited  for  "gross" 
errors.  The  data  is  then  smoothed  and  fitted  to  a  polynomial  and  again 
edited.  The  third  level  of  editing  rejects  data,  either  smoothed  or  un¬ 
smoothed,  that  have  standard  deviations  greater  than  some  specified  criterion. 
The  editing  is  accomplished  in  three  steps. 


■ 
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1)  "Gross"  errors  are  detected  by  comparing  the  change  in  the  measured 
altimetry  (Ah),  between  succession  data  points,  with  the  correspond¬ 
ing  change  in  time  (At).  This  is  a  measure  of  the  altimetry  rate 
(h)  which  is  a  smooth  function.  If  the  altimetry  rate  exceeds 

the  maximum  value  the  point  is  rejected.  The  maximum  for  ft  may 
be  estimated  empirically  from  several  sets  of  GEOS-3  altimetry 
data. 

2)  The  second  level  of  editing  uses  a  n^1  order  polynomial  to  smooth 
the  altimetry  data.  Residuals  are  computed  from  the  difference 
between  the  smoothed  altimetry  and  the  measured.  If  the  residual 
exceeds  10  times  the  standard  deviation  of  the  measurements  the 
point  is  rejected. 

3)  The  third  level  of  editing  rejects  data  based  upon  an  input  sigma 
criterion.  The  criterion  will  be  based  upon  realistic  expected 
accuracy  of  the  GEOS-3  altimeter  measurements.  This  value  will 
be  adjusted  as  experience  is  gained  with  GEOS-3  altimetry  data. 


Plot  Programs 

The  plot  programs  perform  on  option  satellite  ground  track  plots,  land 
outline  plots,  contours,  and  either  geoid  profiles  or  residual  profiles. 

The  land  and  island  areas  are  plotted  from  digitized  land  coordinates. 

The  normal  procedure  for  formal  presentation  of  the  continent  outlines  is 
to  plot  the  ground  tracks  over  a  reproduced  accurate  map. 

The  plot  programs  are  merged  into  a  single  program  called  WMAP.  This 
program  WMAP  consists  of  several  independent  functions  that  may  be  selected 
on  option.  These  options  are  described  below. 

a.  Land  Outlines  -  The  land  outlines  of  the  earth  may  be  plotted  in 
combination  with  either  the  ground  track  or  the  contour  plots. 

b.  Ground  Track  Plot  -  This  program  plots  the  ground  track  of 
satellite  passes  over  the  globe  from  ephemeris  latitude  and  longi¬ 
tude.  Additional  parallels  and  meridians  may  be  selected  with 
program  input  parameters. 

c.  Contour  Plots  -  Contours  of  either  the  standard  deviations  or  the 
actual  geoid  heights  may  be  plotted.  The  standard  deviations  are 
computed  from  the  covariance  matrix  from  the  SARRA  reductions.  The 
geoid  heights  are  computed  from  the  surface  coefficients  saved  from 
SARRA  solution.  In  either  case  the  computations  are  performed  for 
a  grid  of  surface  points.  The  program  includes  the  option  to 
process  unordered  geoid  measurements.  The  original  program  required 
measurements  at  locations  equally  spaced  in  longitude  and  latitude. 
This  later  grid,  in  essence,  outlines  a  rectangular  area.  The 
option  is  of  the  form  that  the  unordered  measurements  are  searched 
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to  provide  points  approximating  those  on  the  grid.  For  most  of 
the  grid  these  points  are  available  and  the  corresponding  contour 
plot  is  obtained.  A  three-dimensional  plot  of  the  geoid  height 
over  the  North  Atlantic  region  was  produced  (see  plotting  attach¬ 
ments)  . 

d.  Residual  Profile  Plot  -  The  altimeter  observational  residuals 

(measured  minus  computed  height)  obtained  from  the  SARRA  reductions 
are  plotted  for  the  purpose  of  visually  reviewing  measurement 
characteristics.  The  residual  profile  plot  provides  immediate 
display  of  altimeter  measurement  noise  and  such  oceanic  details 
as  sea  state  variations.  The  resulting  output  from  this  effort 
was  published^*5)  and  communicated  by  the  Problem  Initiator  to 
other  interested  personnel.  Several  plots  depicting  this  output 
are  attached  for  reference. 
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GROUND  TRACKS  OF  1253  GEOS-3  PASSES 


LONGITUDE 


CONTOURS  OF  GEO  ID  UNDULATIONS  (IOm  INTERVAL)  FROM  THE  GLOBAL  ADJUSTMENT  OF  SATELLITE  ALTIMETRY  AND  GRAVITY  ANOMALIES 


Figure  2.  Contours  of  Geoid  Undulations  (10m  Interval) 
from  the  Global  Adjustments  of  Satellite 
Altimetry  and  Gravity  Anomalies 
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Figure  4.  North  Atlantic  Gravity  Anomalies  from  GEOS-3  Altimetry 
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MATRIX  ELEMENTS  FOR  BOUND/ FREE  THREE  BOW  STATES 
WITH  COULOMB  INTERACTION 


InAjMvton:  Vk.  J.  J a6pe/u,e 
Pfioblm  No:  4932 
PKoje-cX.  No:  1663 

The  L  integrals,  or  matrix  elements,  of  bound/free  three  particle 
states  are  defined  by  Jasperse*-1^  in  his  study  of  three  body  problems.  Com¬ 
putation  of  the  required  volume  integrals  for  the  case  of  a  Coulomb  inter¬ 
action  by  quadrature  has,  for  the  general  case,  proven  unfeasible  due  to 
the  amount  of  time  required  to  achieve  a  useful  degree  of  precision.  The 

work  described  here  utilizes  a  method  of  contour  integration  developed  by 

•  (2) 

L.  Calabi^  ,  combined  with  a  set  of  algorithms  which  allow  the  computer  to 
handle  explicit  high-order  derivative  formulas.  In  this  way  the  inherent 
limitations  of  quadrature  approximation  are  avoided;  leaving,  however, 
truncation  problems  which  may  hopefully  be  adequately  met  by  exploiting  the 
longer  wordlength  capabilities  of  other  computer  systems. 

The  contour  integration  procedure  begins  with  a  breakdown  of  the 
integral  for  a  matrix  element  into  a  sum  of  more  elementary  integrals,  called 
T-integrals,  each  of  which  is  evaluated  by  contour  integration.  One  of 
the  coordinates  of  the  T-integrals,  the  z  coordinate,  is  easily  integrated 
out.  The  remaining  two  integrations  lead  to  "residues  of  residues",  or 
double  residues,  expressed  as  derivative  formulas.  The  computations  described 
here  allowed  the  computer  to  use  explicit  functional  forms  of  the  derivatives, 
carrying  all  quantities  to  the  precision  allowed  by  the  word-length  of  the 
computer.  Nevertheless,  although  there  are  no  numerical  approximations 
involved,  it  was  not  obvious  that  some  truncation  would  not  occur,  there¬ 
fore,  checks  of  both  the  individual  residues  (evaluations  of  T-integrals) 
and  their  sum  (the  L-integrals)  were  devised.  Individual  residues  of  all 
orders  were  also  computed  by  a  symbolic  manipulation  program,  also  utilizing 
full  word-length  precision,  and  agreed  with  those  computed  by  the  above 
algorithms  to  at  least  12  places.  However,  the  range  of  magnitudes  of  the 
residues  was  so  large  that  for  energy  states  higher  than  2  there  was  virtually 
no  quaranteed  precision  on  the  CI)C-(j600  computer.  Since  the  computations 
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involve  complex  numbers,  only  single  precision  is  available  on  the  CDC-6600. 

It  has  been  considered  to  run  the  program  on  an  IBM  machine  with  a  FORTRAN-H 
extended  compiler,  which  is  capable  of  providing  quadruple  wordlength 
precision  with  complex  numbers,  or  approximately  double  the  wordlength  of 
the  CDC-6600. 

All  matrix  elements  for  the  first  and  second  energy  levels  have  been 
computed  and  verified  to  the  precision  of  the  quadrature  method  used  for 
comparison,  and  the  computer  itself.  This  would  seem  to  constitute  a  verifi¬ 
cation  of  the  residual  approach  to  the  problem,  and  also  the  analytical 
details  involved  in  generating  the  computations  to  be  programmed. 

In  his  study  of  applying  contour  integration  methods  to  the  volume 
integrals  developed  by  Jasperse  in  connection  with  the  three  body  problem, 
Calabi's  aim  was  to  devise  a  way  of  computing  these  integrals  with  greater 
precision,  or  at  least  in  less  time,  than  was  possible  with  numerical 
methods.  The  result  of  Calabi's  study  was  a  number  of  derivative  formulas 
which  require  either  special  algorithms  to  evaluate,  or  a  basically  new 
programming  approach  using  symbolic  manipulation  methods.  The  present  effort 
has  led  to  a  package  designed  around  a  set  of  FORTRAN  programmable  algorithms 
for  handling  the  derivative  formulas  in  a  non-numerical  way,  thus  avoiding 
the  limitations  of  numerical  approximation. 

In  addition  the  use  of  symbolic  manipulation  packages  was  explored  as 
a  possible  alternative  to  one  based  upon  FORTRAN  algorithms  and  as  a  means 
of  checking  the  precision  of  the  FORTRAN  computations.  A  symbolic  manipula¬ 
tion  package,  called  SYMBAL,  designed  for  the  CDC-6600  was  eventually  used 
for  checking  purposes,  but  was  too  primitive  in  the  form  used  to  be  integrated 
into  a  useful  package  for  handling  the  L-integrals.  Another  symbolic 
manipulation  program  designed  by  a  member  of  the  Laboratory  was  also  useful 
for  checking  purposes. 

This  report  is  divided  into  two  main  sections.  The  first  section 
presents  a  summary  of  Calabi's  analysis,  which  expanded  the  L-integrals 
into  functions  which  could  be  readily  evaluated  by  contour  integration. 

The  second  section  includes  a  presentation  of  the  algorithms  used  in  the 
computations,  and  the  manner  in  which  the  residues  indicated  in  Calabi's 
report  were  adapted  to  the  algorithms  at  hand. 


Reduction  of  L-Integrals  to  T-Integrals 

Calabi  has  derived  in  detail  the  expansion  of  the  L-integral  into  a 
sum  of  T-integrals,  defining  the  T-integrals  and  the  coefficients  of  the 
expansion,  called  E  coefficients.  The  discussion  in  this  report  will  be 
confined  mainly  to  the  results  of  that  derivation,  omitting  the  develop¬ 
ment  except  where  it  could  help  to  clarify  the  development  of  the  software 
package.  In  several  cases  it  was  necessary  to  re-work  the  derivative 
formulas  to  make  them  compatible  with  the  capabilities  of  the  derivative 
algorithms. 


The  L-integral  is  defined  as: 


-  /‘a.  f\rf 


dz  F(x,y,z)  , 


1  2  2  2  2  2  2 
where  F(x,y,z)  =  x  y  [x  +  2exyz  +  e  y  +  a  ] 


• S  (y)  f  O)  /  ([x2  +  2exyz  +  e2y2]  ^ 

n'"SL'  jt\"!L  Jn'l'  \  ' 

•r([/ 


2  __  ,2  2 
+  2fxyz  +  f  x 


]I/2) 


V 


[x2+2exyz+e':!yZ] 


xz  +  ey 
2  2- 


fx  +  yz 


1/2 


l [y2+2fxyz+f 2x2 ] 


1/2 


Here  the  S's  are  Sturmian  functions  for  the  Coulomb  potential  and  the  P's 
are  Legendre  polynomials.  The  constants  are  all  real  and  verify 

0  _<  Si  <  min  (n,n")  , 

0  H'  <  min  (n'.n'"’)  , 
ef  <  1  , 


l,  all  integers  , 

Following  Calabi 's  notation: 
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where 


N°  ,  z2(p‘«  fa"-»-o;]1/2  „ 

nJl  L  ir(n+&) !  J 


n-i-1 


= 


T  O  n+l 
,2  2.  p=o 

(x  +a  )  r 


E  C(n,Jt,a,p)  x 


2p 


C(n,£,a,p)  =  a 


„  „2(n-*-l-p)  n'£1  r>pl  ^  ,.„prp2^1^  ^  ^  P‘j 


z  c !  2  c_1) 

P1=o  v^nJl  P2=0 


and  is  the  coefficient  of  x*3  in  the  ultraspherical  polynomial 

[see  e.g.,  M.  Abramowitz,  I.  A.  Stegun,  Handbook  of  Mathematical  Functions, 

Nat.  B.  Standards,  Appl .  Math  Series  55,  1964]. 


In  addition,  a  standard  Legendre  polynomial  form  is  adapted  to  the 
arguments  in  the  present  case  as: 

P^)  =  v_A  Z  C(A,p)  uV'P  . 
p=o 

It  is  proved,  without  the  condition  ef<l,  that  the  L-integral  exists 
for  the  domain  of  integration  indicated.  It  is  apparent  from  the  definition 
of  the  Sturmian  functions  that  F(x,y,z)  is  a  product  of  rational  fractions 
of  polynomials.  When  these  functions  are  substituted,  and  the  numerator 
expanded  into  powers  of  products  of  x,  y  and  z  the  L-integral  assumes  the 
form: 


_  oo 


ZEE  E(k,r,s)  T(x,y , z,k,r,s) 
k  r  s 


where 


T(x,y, z,k, r, s)  = 


s  r  k 
x  y  z 


/  2  ,2 ia  (  2  x2\b  /  2  .  2  2  2\C  2  ..  _  2  „2  d 

(x  +6  /  (y  +6  )  [x  +2exyz+e  y  +a  ]  y  +2fxyz+fx  +B 


Also:  a  =  n"+l,  b  =  n""+l,  c  =  n',  d  =  n+l 

N  =  y  NY._  _  N6.,.  N®  .  N3. 

4n  l  n  £  n  1  n£ 


i 


The  function  T  may  be  integrated  over  z  in  closed  form  to  give  a  function 
of  x  and  y  which  is  separated  as  follows: 


T?(x,y)  =  f (x,y)  +  £  fA(f,y)  +  g(x,y)  . 

£  A 

It  will  be  useful  to  present  the  above  functions  in  an  abbreviated  form  to 
show  certain  essential  features  which  determined  the  development  of  the 
software  package.  Departing  somewhat  from  Calabi's  notation: 


f(x,y) 


H  (x,y)  [Ln(y+A)  +  Ln(y-B)  -  Ln(y-A)  -  Ln(y+B)] 

(^2).  (yV)b 


g(x,y)  = 


HQ(x,y)  [Ln(y+A)  +  Ln(y-B)  -  Ln(y-A)  -  Ln(y+B)] 

(x2.«2f  GWf  l wf 


fA(x,y)  = 


(-1)  Hp (x,y) 


(x2*«2)  (yW)  (y2-*2) 


m-r. 


[(y-(-i)Auy)  (y-(-i)Auy)]r^ 


The  functions  A(x)  =  x/e  +  ia/e  and  B(x)  =  fx  +  iB  are  complex  with  the 
bar  indicating  complex  conjugation.  The  function  U^(x)  is  either  A(x)  or 
B(x)  depending  on  the  index  y.  The  indices  X  and  y  are  related  to  r^  which 
runs  from  0  to  d,  with  m=c+d-l.  The  quantities  c  and  d  are  exponents  intro¬ 
duced  above.  The  ((>  appearing  in  the  denominator  factor  in  all  three  functions 
is  itself  a  function  of  x: 


4>(x) 


2  2 
f  2  fa  -eg 

—  x  +  - - — 

e  eg 


2  2 

This  T 2  integral  is  understood  to  be  evaluated  in  the  limit  as  xy(y  -<|>  )  0 

if  necessary  to  guarantee  existence  at  all  points  x  and  y.  Finally  H  (x,y) 
and  H^Cx.y)  are  polynomials  in  x  and  y,  the  orders  of  which  depend  on  the 
exponents  k,  r  and  s. 


Denoting  the  Cauchy  principal  value  by  CPV,  it  is  proved  that: 


CPV 


dy  T2(x,y) 
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2  2 

exists  if  fa  >eB  ,  and  may  be  evaluated  as  the  sum  of  the  residues  of  f,  g 

and  fX  as  long  as  the  condition  max  (n',n+l)  <_  Z  +  l'  +  5  is  satisfied. 

The  latter  condition  is  the  only  constraint  on  the  basic  parameters  of  the 

2  2 

problem.  The  relationship  fa  <eg  may  be  accommodated  by  changing  the  order 
of  integration  or  interchanging  variables. 

The  above  integral  is  carried  out  as  a  succession  of  two  contour  inte¬ 
grations.  First  the  quantity 

T3(x)  =  CPV  J  T 2(x,y)  dy 

-CO 

is  computed  as : 

T3(x)  =  27riZ  Res  [T2(x,y.)]  . 

i 

Since  T2(x,y)  is  broken  into  several  components,  with  contours  chosen 
independently  for  the  several  components,  the  y^  are  whatever  poles  are 
necessary  to  include  in  all  of  the  contours  Likewise, 


~oo 


=  2iri£  Res  [T^fx.)]  , 

j  '  J 

or 

T  =  -4tt2  Z  I  Res  [Res[T2(x  ,y.)]]  , 

j  i  J 

leading  to  a  double  residue  computation. 

The  double  residue  in  the  expression  for  T^  is  equivalent  to  a  double 
derivative  in  x  and  y,  of  each  of  the  component  functions  of  T2(x,y).  The 
order  of  the  derivative  in  each  case  is  determined  by  the  order  of  the  pole 
in  each  variable.  It  is  immediately  apparent  from  the  functions  f,  g  and 
fX,  and  from  the  potentially  great  number  of  T-integrals  in  an  L-integral 
expansion,  that  computer  evaluation  of  the  residue  derivatives  is  necessary. 

It  also  becomes  clear  from  the  order  of  the  derivatives  and  from  the  precision 
required  that  any  type  of  numerical  differentiation  would  be  inadequate. 
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Calabi  develops  formulas  to  evaluate  the  required  derivatives  as 
programmable  algebraic  expressions  which  avoid  the  precision  limitations  of 
numerical  differentiation.  These  expressions  are,  however,  extremely 
cumbersome  and  were  replaced  by  derivative  algorithms  which  are  basically  a 
much  more  compact  way  of  expressing  the  above  algebraic  evaluation.  In  most 


cases  using  the  algorithms,  it  is  a  relatively  simple  matter  to  code  an 
individual  residue  computation,  although  in  some  cases,  to  be  seen  below, 
considerable  expansion  of  the  derivatives  is  necessary  before  the  algorithms 
can  be  applied. 


Calabi's  algebraic  expressions  also  permit  an  analysis  of  the  existence 
of  individual  T-integrals.  While  the  existence  of  the  Cauchy  principal 
value  of  the  L-integral  is  proved  without  restriction,  the  existence  of  all 
T-integrals  without  condition  is  not  guaranteed.  This  is  due  to  the  fact 
that  the  expansion  of  the  L-integrals  can  introduce  special  features  into 
the  component  T-integrals.  In  particular,  Calabi  shows  that  the  contour 
integral  leading  to  one  of  the  residues  exists  only  under  the  condition  noted 
above  that  max  (n%n+l)  <_  l  +  l'  +  5.  (See  pp.  32r-33).  Whether  this 
condition  could  be  removed  or  relaxed  by  a  suitable  manipulation  of  the 
functions  may  be  a  matter  for  further  study. 

Although  the  individual  poles  in  the  T2(x,y)  function  are  not  difficult 
to  discover  and  isolate,  being  the  roots  of  explicit  linear  and  quadratic 
forms,  their  distribution  presents  some  difficulty  in  view  of  the  branch 
points  of  the  log  functions.  Calabi  develops  an  integration  theorem  for 
the  complex  plane  which  gives  a  procedure  for  grouping  the  components  of  the 
T2(x,y)  function  into  different  regions  of  the  complex  plane  in  accordance 
with  their  poles  and  branch  points. 


Looking  first  at  the  y  integration,  the  resolution  of  T2  into  f,  g,  and 
fA  provides  that,  for  x  real,  the  contour  integrals  of  f  and  fA  may  be 
taken  in  the  positive  imaginary  part  of  the  y+iv  plane,  with  the  branch 
points  restricted  to  the  negative  half  of  the  plane  (including  v=o) .  Similarly, 
g  may  be  integrated  over  the  negative  y+iv  plane,  with  branch  points  restricted 
to  the  positive  half  of  the  plane. 


It  may  be  observed,  for  example,  that  in 


f(x,y)  = 


HQ(x,y)  [Ln(y+A)  +  Ln(y-BJ  -  Ln(y-A)  -  L(y+B)] 

(y wf  (y2-t2)" 


the  log  functions  have  branch  points  only  in  the  negative  imaginary  part  of 
the  y+iv  plane,  for  real  x,  since  the  imaginary  parts  of  A  and  B  are  positive. 
Further,  since  by  the  integration  theorem  it  is  only  necessary  to  include 
poles  in  the  v>0  part  of  the  plane,  excluding  the  real  axis,  the  only  pole 
in  the  y  plane  that  must  be  included  is  that  at  y=+iy,  of  order  b.  In  this 
way: 


CPV 


J  f(x,y)  dy  =  2iri  Res  [f(x,iy)]  , 


where 


Res  [f (x, iy) ]  = 


(b-l)l 


-b-1,  ..b 

3  (y+y) 


3y 


b-1 


f(x,y) 


y=iy 


If  the  same  integration  theorem  is  invoked  for  a  contour  integration 
in  the  x+iw  plane,  it  becomes  apparent  on  inspection  of  the  arguments  of 
the  log  terms  that  two  of  the  log  terms  have  branch  points  in  the  upper  half 
of  the  x+iw  plane  [i.e.,  iy-B  =  iy  -  fx  +  iB,  and  iy  -  A  =  iy  -  —  +  i  — ] , 
and  the  other  two  have  branch  points  in  the  lower  half  of  the  x+iw  plane. 

If  the  function  f  is  further  broken  down  to  isolate  these  branch  points,  with 
f^(x,y)  and  f2(x,y)  having  branch  points  in  the  lower  half  and  upper  half, 
respectively  of  the  x+iw  plane,  then  the  above  residue  computation  becomes: 


CPV 


f (x,y)  dy  =  2iri 


Res  [f L (x, iy) ]  +  Res  [f2(x,iy)] 


Avoiding  those  parts  of  the  x+iw  plane  with  branch  points  it  is  apparent 
that  Re^f^(x,iy)]  must  be  integrated  in  the  upper  half  of  the  complex  plane 
and  Res[f2(x,y)]  in  the  lower  half.  It  is  also  clear  that  both  Resff^]  and 
Res[f2]  have  poles  at  x  =  +i6  of  order  a.  From  the  above  definition  of 
<j>(x)  it  is  also  clear  that  upon  substituting  y=iy  after  differentiation,  the 
denominator  factor  produces  the  quantity  ( - y^ ^ jn  the 
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■  >-■  ■■*  . 


2  2  f 

denominator.  Then  by  factoring:  4>  (x)  +  y  =  ^-(x+iv)  (x-iv) ,  where  v  is 
a  positive  constant.  Thus  both  f^  and  also  have  poles  in  the  x  plane  at 
x  =  ±iv.  Thus  the  contour  integrals  of  Res  [f ^ (x, iy) ]  and  Res [f 2 (x, iy) ]  must 
include  residues  at  x  =  i£,iv  and  x  =  -i6,-iv  respectively.  Thus: 

/oo  r  °o 

J  f(x,y)  dxdy  =  -4tt2  ^  Res  [Res  [f  ^  (i<$ ,  iy)  ]  ] 

— cx>  _oo  ' 

+  Res [Res [ f ^ (iv, iy)  ]  ] 

+  Res[Res [f2(-i6,iy)]  ] 

+  Res(Res[f2(-iv,iy)]]  |  . 


Res [Res [fj (i6,iy)]]  =  Cb_i) ! (a-l) 


a-1  b-1 

~ — r  — bT  ^x"i6^  Cy-i-vD  f,(x,y) 
Sx3'-1  8y  1  1 


.  am+b-2  , 

Res  [Res  [f  x  Civ ,  iy)  ]  ]  =  ^ ,  (ii»b-2) !  7m^2  (x-iv^+  ' 

dX 


— g-(y-iy)  f^x.y) 


In  the  expression  for  Res [Res [f ^ (i£ , i y) ] ] ,  where  the  poles  in  both  the 
y  plane  and  the  x  plane  may  be  factored  out  of  the  denominator  before  any 
derivatives  are  taken,  the  algorithms  make  it  possible  to  code  the  double 
derivative  by  simply  calling  general  subroutines  as  described  below.  In 
the  expression  for  Res [Res [f  (iv, iy) ] ]  the  x  pole  may  not  be  factored  out 
of  the  denominator  until  the  y  derivative  has  been  taken,  so  that  some  pre¬ 
liminary  manipulation  must  take  place  before  applying  the  double  derivative 
algorithm.  This  is  accomplished  by  writing 


fjU.y)  = 


_ P(x,y) 

(y2+y2)  (y2-4>2) 


2  2  ~m 

and  by  expressing  the  derivative  and  extracting  the  factors  of  (y  -<(>  ) 
in  a  general  form  before  the  derivative  algorithms  are  applied. 
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Differential  Algorithm  Formulations  for  the  Residue  Evaluations 

The  values  for  the  residues  necessary  to  determine  the  atomic  Coulomb 
integrals  are  determined  by  incorporating  algorithms  for  evaluating 
derivatives  of  functions.  This  procedure  is  equivalent  to  a  straight  forward 
computer  evaluation  of  mathematical  formulations,  thereby,  having  accuracy 
proportional  to  the  computer  word  length.  A  procedure  of  this  form  is  a 
necessity  for  this  problem,  in  that,  approximation  by  finite  difference 
methods  is  too  inaccurate.  The  two  basic  algorithms  used  are  the  following: 


Let  the  derivatives  with  respect  to  x,y  exist  for  functions  F(x,y) 
and  G(x,y) ,  then 


Derivative  of  a  product  of  functions 
DI,J  (F(x,y)G(x,y))  =  Z  Z 


x..y 


1=0  j=o 

Derivative  of  a  reciprocal  of  a  function 


F(x,y)  D^’J-jG(x,y)  .  (Al) 

J/  A>7  A>7 


D'l^  (1/F(x,y) ) 

A  >  7 


(A2) 


In  regards  to  the  computational  needs  for  this  problem,  the  above 
algorithms  are  the  fundamentals  from  which  other  needed  algorithms  were 
constructed.  For  example,  the  later  algorithm  is  needed  for  constructing 
a  formulation  for  evaluating  derivatives  of  the  Natural  Logarithm  of  a 
function  F(x,y),  that  is, 


D l’*  (An(F (x,y))) 
A  y  y 


I 

z 

i=o 


D1  lTDI-i-J-j 

x,y  F(x,y)  x,y 


F(x,y) 


(A3) 


Note:  J>o;  for  J=o 


D°’°  =  «-nCpCx,y))  and 

A>  7 


DI,Q  («,n(F(x,y)))  =  Z  f1"1)  D1,P  rr-  1  .  D1'1*0  F(x,y) 
x,y  ,7JJJ  i=Q  \  1  /  x,y  F (x,y)  x,y  v  ,7J  ‘ 


Another  formulation  required,  of  a  fundamental  nature  for  this  problem,  is 
to  evaluate 
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Di:Jg  p(x'y) 


y=4>(x) 

x=x 


2  1/2 

where  4>(x)  =  (Ax  +  B) 

2  2 

and  p(x,y)  =  PQ  +  Pj  y  +  P2y  +  P3xy  +  p4x  +  p£x 

a  quadratic  polynomial  in  x,y.  It  was  decided  to  treat  this  as  follows: 

Take  the  y  differentials  and  substitute  y=<t>(x),  thereby  requiring  the 
evaluation  of 


Dx,y  p(x’y) 


y=<l>00 

x=x 


D*  p(x,<Kx)) 


x=x 


D*  l  p(x,y) 


y=<p(x) 

x=x 


(p1  +  2p2<Kx)) 


x=x 


Dx,y 


y=4>(x) 

x=x 


D*  (2p2) 


x=x 


Dl,Jy  P(X-rt 


o  for  J>2 


y=<t>(x) 

x=x 


Then  algebraic  differentials  of  the  irrational  function  <j>(x)  were  obtained 
manually,  for  example: 

1/2 


D°  <t>(x) 


(ax2+b) 


Dx  c(.(x)  Ax 


(ax2+b) 


-1/2 


2  /  2  \ ~^/^  22  2  -3/2 

<J>  (x)  -*■  A^Ax  +B)  -  A  x  (Ax  +B) 


etc. 


and  thus  the  required  formulation  was  accomplished. 


(A4) 
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The  first  (2)  of  each  of  these  formulations  is  utilized  in  each  residue 
calculation,  whereas,  the  later  (2)  is  required  only  for  certain  of  these 
residues.  What  follows  is  in  general,  a  differential  formulation  for 
specified  residue  evaluations. 

Residue  (x=i6,y=iY)  (including  fA) 

These  are  the  simplest  of  all  the  required  residue  calculations. 

They  are  obtained  by  a  direct  application  of  equations  A1 ,  A2,  and  A3. 


Residue  (x=iv,y=iy)  (not  including  fA) 

These  and  subsequent  residues  are  of  a  nature  that  upon  performing 
the  y  differential  corresponding  to  the  order  of  that  pole,  the  order 
of  the  x  pole  is  thereby  changed.  Hence,  in  addition  to  utilizing 
equations  A1,A2,  and  A3  the  following  formulation  must  be  implemented 
into  the  calculation  and  upon  taking  the  y  differentiation  the  resulting 
equations  to  be  programmed  are: 


n  =  1  , 


3X  (x+iv) 


- — j-  [(x+iv)  (x-iv)  Py-m(2y)  P] 


y=iy 

x=iv 


=  2 


1 


^  in  ,  . 

ax  (x+iv) 


i72  CCX+iV)2  (X-iV)2  Pyy 


-2m(x+iv)  (x-iv)  (2y  Py+P) 


+  (m+1)  (m)  (2yr  P] 


y=iY 

x=iv 


where  P  is  the  functional  form  of  the  "reduced  integrand"  with  the  sub¬ 
scripting  P  ,  etc.,  representing  differentiation  with  respect  to  y. 
Currently,  this  is  the  formulation  used  in  the  program;  however,  a 
general  form  for  n  is  available  and  has  been  separately  programmed, 
which  is : 
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Residue  (iy.iv)  = 


tWQ 


(b-1) 


b-1  i  i-y  r 
III  l 

1=0  Y  q2 


(m+y-1) ! (m+i-y-1) !  ,  .  H2 


(m-1) !' 


(-D 


,  .  ,1_qrq2 

(-iy) _ nm 

(m+i-1)!  x 


ql+q2 

Dm+i-l  D  b-1  -i  ((>  P(x,y) 

^  (-f/e(x+ispv) )m+* 


|x=ispv 

y=iy 


Residue  (x=i6,y=±<))) 


In  brief,  the  computation  of  the  <j>  residues  requires  the  simultaneous 
consideration  of  g(x,y)  at  y  =  ±<J>  in  the  form 

G(x)  =  Res(x,<t>)  +  Res(x,-<{>) 

Care  must  be  taken  in  choosing  the  planes  of  integration  to  avoid  conflict 
with  the  branch  points  from  the  logarithm  term.  Upon  substituting  y=<t>, 
the  Ln  factor  has  the  form 

2  2 

Ln  u  =  Ln  [a  -y  ]  =  Ln  [E(x+in+)  (x+in  )] 
or 

Ln  [E(x+in+)]  +  Ln  [x+in  ] 

wherein  the  sign  of  n  ,  o  determines  the  contour  of  integration.  Ex¬ 
pressing  the  integral  g  as  the  sum  of  integrals  over  G^(n+)  and  G2(n  ) 
the  following  formulations  were  determined  for  these  residues: 

1  a-1  m-1  Ln  tE(x+iTV)] 

Res(G. ,i6-sp)  =  7 — prf-p — r-r~r-  D  D  - 

1  (m-1)! (a-1)!  x  y  ,  .  .  .a 

J  J  (x+i<$*sp) 


P(x,y) 


y  =  $ 

x  =  iS*sp 


Res(G  ,i6*sm)  =  7 — pyj-p — ppp-  D  3-1  D  m_1  — LMuJ - 

2>  (m-1)! (a-1) !  x  y  (x+i6.sm)a 


P(x.y) 


x  =  i6*sm 
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where  u  is  loaded  at  D°  by  [x+in  ]  and  sp  =  sign  (n+) ,  sm  =  sign  (n_) » 
and  E  =  constant.  Hence,  these  residues  were  evaluated  by  programming 
this  formulation  in  conjunction  with  equations  A1 ,  A2,  A3,  and  A4. 

Below,  just  the  Residue  Formulations  are  formally  outlined  with 
partial  mention  of  referencing  to  equations  Al,  A2,  A3,  A4. 

Residue  (x=iv,y-±<t>) 

The  mathematical  formulation  derived  for  this  residue  is 


m-1 


Res(G1,iV  sp)  -  (m+b-2) !  ^  (  k  5  ^m-k-P 


(— )  (m-k-1)!  D  b+k_1  (x+ivsp)  b_k 

6  X 


•  Ln  [E(x+in+)]  Dym'1_k  P(x,y) 


k 

•  E  W(b,j,k) 
j=o 


(<t>+iY)k_^ 


(•My)’’ 


y 

x 


4> 

iv  *sp 


W(b, j ,k) 


(-l)k  (k) 


(b+j-l)l  (B+k-j-1)! 
[ (b— 1 )  !  ]2 


Residue  (x=-in+ , y=±<J> ) 

This  residue,  precipitating  from  the  differentiation  of  the 
logarithmic  term,  is  as  follows: 


Res(G2,-in+)  = 


,  m-1  £  .  „ 

_ i _  E  g  fm-^)  f*") 

Cm-2)!(m-l)!  £=1  j=Q  (  »  3  V 


C-l-*.)  ("»-!-&)  !  («.-!) !  (-1)A_1  Cl+(-l)£‘j)  E"£ 


^£-1  ,  .  . -£  .j  £-j  nm-£-l 
D  (x+in  )  <t>J  a  D  H(x,y) 

x  "  y 


y  =  <p 
x  =  -in. 
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r 


La 


Residue  (x=in,y=<|>=o) 

This  residue,  upon  considering  an  x  pole  corresponding  to  <J>=o ,  has 
the  following  formulation: 


Res(<j>=o)  =  const 


»->  s  £  W; 


2"m"Jl+j  D™"1  Ln(E(X+ip+) ) 
(IX- 1)1  (E(x+iwsp))!) 


1  (-l)*"j (r-k-l)i 


(r-k-l-j) ! (m-1) ! 


.r-k-l-m-£  ..m-l-d  ,  . 

(>  Dy  P(x,y) 


y=4> 

x=iw, (n+,n-) 


(A-i)l 


n  IX-l  l-£ 

D  a, 

x  L 


(IX-l)!  2  2, 

aL  7  (2-(x+iwsm)) 


,r-k-l-m-q+e  ,  „ 

♦ - 

IX  y  ^q,y 


y  =4> 

x=iw(n) 


Dy(Px,y) 


J(S-k_1 

(x2+62) 


-k-1  k  min(t , 2 j )  t^s  n  v  min (2(k-j ) ,t-s-u) 

7^  £  L  Z  E  E  zl 

j=o  s=j  n=o  y=o  q-=o  r=o 


Sn^v 

£ 

D  qi =072,4,  < 


s=J 

t-s-n-r 


n=o  y=o  q^o 


,  -J  f  2  ^-j  /k 
(eg)  (e  ) 


>2=o  qi=U,  2,4,6  P^OTM.h 

or  1 ,3,5,7  or  1 ,3,5,7 


0(77) 
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(;)  (:)  (:)  (:;)  (;)  (;)  (t) 


(:,) 


(j  !)  (b+v-1) !  (b+n-v-1) !  (k-j) ! 

(2j -s) !  (b-1)!2  (k-j -r)  :  (k-j -t+s+n+r) ! 


22j  -S(-1) 


(q2+p2  +  (q2+q!+p2+pp/2) 


(tT  („T  (>WY°  (*W) 


2  2\k_^  -t+s+n 


t+2j -2s-2p-2q 


W  =  V"Fi-  -  T  q!  +  q2  =  2q 


n2  =  (x2+a23/e“  p1  +  p2  =  2p  . 

It  should  be  noted  that  the  residues  described  above  are  the  so- 
called  f  and  g  residues  or  non-X. 

The  following  section  describes  the  X  residues  and  it  should  be 
noted  that  of  the  four  basic  formulations  described  above  only  A1  and 
A2  are  required. 

The  X  residues  of  form  (x=iS,y=iX)  and  (x=iv,y=iy)  will  not  be 
discussed  since  the  differential  formulations  are  similar  to  those  of  the  same 
notation  described  earlier.  The  main  difference  being  only  in  the 
"reduced  integrand",  i.e.,  what  is  to  be  differentiated. 

Residue  (x-xq(u,uX) ,y=iy) 

The  formulation  for  this  residue  being 


case  b=2 


1  r 

Residue  =  pyyp  ^"(Cx-x^  Py  -  rn  P) 
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case  b=3 


Residue  = 


2(r  +1)!  x 
n 


r  +1  ? 

D  (Cx-x  )  P  =2r  (x-x  )P  +r  (r  JP) 
o  yy  o  y  n  n+1'  1 


y=iy 

x=x 

o 


Residue  (x-xQ(y,Y,A) ,y=u,u) 


Basically,  these  residues  are  evaluated  as  follows:  for  n=2,  rn=l,2. 

case  r  =2 
n 


Residue  = 


J. 

b! 


DxC*-xo)Py-p) 


y=u,u 

x=x 

o 


Residue  (x=xo(y,<J>,A)  ,y=u,u) 


case  r  =1  : 
n 


Residue  =  Dm~2 

(m-2) !  x 


y=n,u 

x=x 

o 


case  r  =2: 
n 


Residue  = 


1 

(m-2)! 


D™~2((y2-<f>2)p-(m-2)(2y)P) 


y=n,u 

x-x 

o 


As  can  be  expected  in  a  problem  of  this  magnitude,  a  number  of  analyti¬ 
cal  requirements  were  incurred.  They  included: 

Time  Estimate  -  In  the  early  stages  of  this  problem,  an  effort  was 
made  to  forecast  the  computer  time  requirements.  Basically,  this  was  done 
by  comparing  the  time  for  a  calculation  of  a  lower  order  integral  and  that 
of  a  higher  order  integral.  This  data  was  then  correlated  to  the  calculation 
procedure,  namely,  the  number  of  residue  calculations  required  and  the 
order  of  each  differentiation.  From  this  some  time  estimating  has  been 
done  but  is  not  assumed  accurate. 
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Verification  of  Residue  Calculation  -  In  addition  to  the  bulk  of 
checking  mentioned  in  the  introductory  section,  it  should  be  noted  that 
some  lower  order  residue  calculations  were  checked  by  hand  and  by  numerical 
approximation. 

Error  Estimate  -  Also  of  primary  concern  is  the  number  of  digits  of 
accuracy  in  an  internal  calculation.  As  noted,  the  residue  evaluations 
(for  n=2),  are  not  the  problem  in  theirself  but,  rather  in  the  summing  of 
them  to  form  the  value  for  the  integral.  In  essence,  the  loss  of  accuracy 
is  caused  by  combining  terms  of  equivalent  magnitude  but  opposite  in  sign. 
This  later  case  can  apparently  be  determined  by  comparing  the  separate 
sum  for  both  the  positive  or  negative  residues.  A  more  general  form  is 
naturally  more  desirable.  Presently,  being  considered  is  to  calculate  the 
same  integral  using  different  word  lengths,  which  when  compared  then  exhibits 
the  truncation  error. 

Results  -  Table  I  contains  a  family  of  Coloumb  integrals 

of  the  form  L  *  for  two  earlier  investigations  and  the 

n,n  ,n  ,n  ,1,1  & 

current  effort  for  data  a  =  2.31,  3=  1.14,  r=6=  2.54555,  E  =  .999864, 

F  =  .05  where,  in  brief,  Jasperse  set  Z  =  l'  =  o  and  then  parallels  Calabi; 
whereas,  Grossbard  utilized  a  strightforward  3-dimensional  numerical 
integration. 

This  effort  for  higher  orders  is  presently  being  continued  as  a  part 
of  another  Air  Force  project. 
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TABLE  I 


Grossbard 

Jasperse 

Calabi 

L111100 

3.86290148436 

3.86265601596 

3.862656015929 

L1 12100 

-.493275047545 

-.493340150938 

-.4933401497509 

L111200 

-3.9007 

-3.90067764136 

-3.900677740736 

L112200 

3.3943 

3.39404022841 

3.3940402366 

L121100 

-2.16382206187 

-2.16297156832 

-2.162971562578 

L121200 

3.9085 

3.90817439123 

3.908174487 

L122100 

4.0236 

4.02337484402 

4.023373012508 

L122200 

-6.1341 

-6.12737514331 

-6.1273690616 

L211100 

.0128271740761 

.0128830519028 

.01288304718 

L211200 

2.6816 

2.68125017719 

2.68125013708 

L212100 

.928635340444 

.928580637072 

.92858049261 

L212200 

-2.2927 

-2.29275562759 

-2.29275790467 

L221100 

3.58132361885 

3.58092539223 

3.5809311925 

L221200 

-5.6336 

-5.63342551987 

-5.6334222233 

L222100 

-1.4738 

-1.47375365940 

-1.4737602310 

L222000 

5.6309 

5.63090879091 

5.630200201 

L121201 

-.322285 

-.32213486 

L1222 01 

-4.06004 

-4.0661064 

L212110 

.767619 

.76769742 

L212210 

-3.58466 

3.5845788 

L221201 

4.14230 

4.1416992 

L222110 

2.8789 

2.8706734 

L222201 

-.679022 

-.67904956 

L222210 

-1.03949 

-1.0390567 

L2222 1 1 

-.0653925 

-.065392053 

i 
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ELECTRON  A NV  PARTICLE  VISTRIBUTION  FUNCTIONS 


Initiaton:  Vn.  J.  JctspeAie 

Problem  No:  40 37  [fiomeAly  Pnoblm  No:  4916) 

Project  No:  464 3 

Over  the  last  three  years  the  mathematical  and  numerical  analyses  has 
been  developed  and  the  computer  programs  written  to  calculate  the  electron 
and  particle  distribution  functions  in  the  ionosphere.  This  was  accomplished 
using  a  non-linear  difference  equation  approximation  of  an  integral  - 
differential  equation  which  described  the  energy  transport  equation  (see 
references  1,2, 3, 4, 5).  Computer  programs  have  been  developed  which  use 
various  approximations  to  a  larger  more  complete  model.  The  complete  model 
has  been  programmed  and  the  results  presented  in  references  4,5.  The  complete 
model  is  known  as  the  BFPEDF  code.  This  model  and  the  numerical  approxima¬ 
tions  developed  by  Boston  College  are  described  below. 


Simulation  of  the  BFPEDF  Code 


BFPEDF  solves  the  following  differential  equation 

/  2me  \  1  ,  . 

^m.  +m  )  Yeis  Z 


y  (z,E)  E3/2  +  Z 
'mn''  1  ^  ..  ...  , 

s  '  is  e' 


r  (^i  \ 

.  |  m  ,+m  I 
J  \  "J  e/ 

f  EM  f  Ew  i  j 

+  K  I  (z,E)  H(z,E)  =  +  /  dE’  Z  G  ..(z,E’)  +  dE1  Z  G  ..(z,E  ) 

ee  o  '  j  pij v  j  enij  v  ’ 

E 

E  i/o  I  -l  m  ■*1/2  _ 

-/  dE1  ECE1)  ’'rs(z-E  )•  H<z-E  ’  -  J  dE  EE<E  >  TtkU.Eh 
E  s 

fl<E‘-Etk> 

-/E"  dE1 


H(z,EA) 

+  J  dE 

E 

ZZ(E  +E  . 
tk  tk 

ZZ  exp  I 
tk  L 

Etk 

}  (E'*Etk 

;  (KTTtk(z)). 

ZZ  exp  I 

tk  *- 

’  Etk 

I  ,  1/2 

J  (E1) 

■  ‘Vtk*1” 

E  tk 

1/2 
I 

1/2 
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(2m  \ 

- “ — ] 

m  .  +m  I 

nj  e/ 

(’  2m  \ 

- —  1  1 

m.  +m  ) 
is  ej 


W‘-E)  3/2  W’> 


Yeis(z)  W*0  +  [VZ’E)  +  J^Z'E^  3FH(Z>E) 


for  H(z,E)  versus  E  for  a  particular  z  will  be  constructed  as  follows: 

TERM1  +  TERM2  +  TERM3  =  TERM4  +  TERM  5 
+  TERM6  +  TERM7  +  TERM8  +  TERM9  +  TERM10 
+  TERM11  +  TERM12  +  TERM13  +  TERM14 

where  the  terms  are  numbered  according  to  how  they  appear  in  the  formula. 
Thus  TERM13  = 

T*  *2U.b)  sfrH-  ■ 

In  addition,  the  following  boundary  condition  was  used  in  the  solution: 


■/' 


M  dE(E1/2)  H(z,E)/ne(z)  . 


The  method  used  to  solve  this  differential  equation  was  to  guess  a  point 

distribution  for  E  s.t.  0  =  E,  <  E0  <  E,  <  ...  <  EMC.,.T  =  E., 

123  NEVAL  M 

For  E  =  E2,E3,...,  eneval  the  equation  was  then  approximated  by  a  non¬ 
linear  difference  equation  with  values  of  H(z,E)  defined  at  E^ ,E2 ,E3> . . . , 
ENEVAL'  a£^ition,  the  boundary  condition  was  also  approximated  by  a 

non-linear  difference  equation. 

These  equations  are  numbered  as  follows: 

1.  boundary  condition 

2.  Equation  for  E2 

3.  Equation  for  E 


NEVAL  Equation  for  E 


NEVAL 


The  unknowns 

in  these 

equations 

ENEVAL' 

Thus  we  wish 

to  solve 

a  matrix  < 

V 

Al,2 

Al,3, 

A2,l 

A2,2 

A2,3, 

J1 ’2’  3’ 


Neval,!,  Neval,  2,  'Sjeval.s, 


....,  A 
A 


1,  NEVAL 

2,  NEVAL 


■  ■  *  A, 


NEVAL, NEVAL 


Hfz.Ej) 

H(z,E2) 


H(z,E. 


NEVAL 


where:  A^  ^  is  the  value  associated  with  H(z,E^.)  for  equation  i.  H(z,E^) 
is  the  desired  solution  and  are  constants.  The  difficulty  in  solving 
this  equation  is  that  the  A^  ^'s  are  functions  of  the  H(z,Ep's  which  makes 
the  equation  non-linear. 


The  method  of  solution  for  the  matrix  equation  is  as  follows: 

Guess  a  value  for  HCz.Ep,  H(z,E2) , . . .  ,H(z,ENEVAL)  say  Hold(z,E1), 
H°ld(z,E2) , . . . ,H°^d(z,ENEVAL)  then  calculate 


ERR°ld 

Al.l 

Al,2 

’  Al, NEVAL 

ERk°ld 

= 

A2,l 

A2,2 

A2, NEVAL 

errSe?al 

aneval,i 

ANEVAL,2, 

•••  ANEVAL, NEVAL 

_  — 

■■  ■  ■ 

- - 

H°ld(z,E1) 

Bi 

H°ld(z,E2) 

B2 

.  .0 Id  ,  P 
,  1  ’NEVAL) 

bneval 

■- 


where  A.  .  is  calculated  using  the  values  of  H  .If 

i.J 


1 

NEVAL 


(ERR°ld)  <  RMSMAX 


where  RMSMAX  is  a  constant,  then  the  H  vector  is  the  desired  solution, 

n  ew 

otherwise,  a  new  H  vector  H  is  calculated  as  follows: 


Cl,l  Cl,2 

Cl,3, 

Cl, NEVAL 

C2,l  C2,2 

C2,3, 

C2, NEVAL 

= 

CNEVAL,1  CNEVAL,2 

CNEVAL,3, • 

‘  -  ’  CNEVAL, NEVAL 

9  ERRX 

9  ERR 

+ 

9  ERRj 

9  H°ld(z,E1)  9 

H°ld(Z,E2) 

’  9  H°ld(z, 

ENEVAL) 

9  ERR2 

9  ERR2 

9  ERRX 

9  H°ld(z,E1)  9 

..old,  „  .. 

H  (z,E2) 

9  H°ld(z, 

eneval^ 

9  errneval 


9  ERR, 


NEVAL 


9  H°ld(z,E1)  9  H0ld(z,E2) 


9  er1Wal 

5  H°ld<z-ENEVAL> 


Then  solve  the  following  matrix  eouation 


'1,  NEVAL 


:1,1 

Cl,2 

2 , 1 

C2,2 

'NEVAL,  1 

CNEVAL, 2 , . . . , 

J2 , NEVAL 


^NEVAL 


ERROld 


erCal 


'  :  •  •; 


and  the  new  guess  for  the  H  vector  is 


\  v 

u 


new 

H(z,El) 

old 

(z.Ej) 

A1 

new 

Cz,e2) 

oid 

(z,E2) 

1 

+  - 

A2 

• 

= 

2* 

• 

Hnew 

^z,enevaiP 

H°ld 

^ z ’ enevaiP 

aneval 

SL  =  0  or  1 ,  or 

2 ... is 

defined  as  the  smallest 

integer  s. 

,  p'JEVAL  2  ,  /NEVAL  ,,  2 

NEW  J  <ERRrW)  <  NEVAL  /  J j  (ERR^  > 


After  Hnew  is  found  H°^d  is  redefined  as  HneW  and  the  method  is  iterated 


To  start  the  iterative  process  we  need  an  initial  guess  for  HCz.E^) 


This  initial  guess  is  usually  provided  by  guessing  Te(z)  the  n^'s  and 


approximating  the  H(z,Ei)  by  Hmax(z,E.)  +  H^z^),  where 


3/2 


W2’E1>  *  2,meC2)  l,ye(z)l  [- 


and 


e)  =  (<4)  s  [‘  ■ e,ip  (■  w*>)J  (h) 


^CSD  ^ 2  ’ E’  ^ 


tk  L1  "  GXP  V.  KrTtk< 


p  )  Ytk(z,E.) 


. 1/2  ZZ  n  •  (z)  0(E.  -E.  .  „)  Q.  , 

1  jfc  \  Ei  /  nj  1  1J* 


eCE.-E^..)  +  K  E 
v  l  tk'  ev 


/  2  m  \  .  3/2  /  2  m  \ 

Z  { - - — \  y  (z,E.)  +  Z  (-4-)  ( - — )  y  •  (z) 

.  Im  .  +  m  I  'mnj  ^  iJ  vl"*  Im.  +  m  /  'eisv  ’ 

J  V.  nJ  e/  s  E  V.is  e/ 


fcf  ■: 


ee 


(z) 


/  dE  Z  G  .  .  (z.E1) 

E  4  PM 

i  J 


-1 


X 


Where  Ytk('Z,Ei')  =  nnt ^  KeV  ^Ep  7  Qtk ^Ei^  and  for  t=i»K=l,7: 


(  0  for  Ej  -  E  <  Etk 

W  =V  A  E  ♦*  E  B-  N" 

q  A^.i  E  .  E  . 

no  tk  (  tk.  M  ,  tk.  c  ^  ^  c 

—7-  (“E7>  -  (~e7 >  '  for  Etk  i  E  i  Em 

vEtk 


with  Tik(z)  =  Tet(t) 


b.  for  t=2,  k=l-8:  Q2j<. C E)  is  read  in  as  a  table  of  C^CE)  versus  E 

values,  where  additional  values  of  Q^^CE),  not  directly  found  in 
the  table,  are  approximated  by  linear  interpolation.  In  this  case 
T2k^  =  TVIB^  E°r  t=^’  ^=9-23:  f°un<i  similarly  to 

the  case  t=l,k=l,7.  In  this  case  T2k(z)  =  TeX,^Z^  ' 

c.  for  t=3,  k=l-8:  Q^^(E)  is  read  in  as  a  table  of  Q^CE)  versus  E 

values,  where  additional  values  of  Q^CE) ,  not  directly  found  in 
the  table,  are  approximated  by  linear  interpolation.  In  this 
case  T3k(z)  =  TVIBU). 


d.  for  t=3,  k=9-18:  ds  f°und  similarly  to  the  case  t=l, 

k=l,7.  In  this  case  T^^Cz)  =  T  . 


e.  for  t=4 ,  k=l,2  and  t=5,  k=l:  Qtk(E)  is  read  in  as  a  table  of 
Qtk(E)  versus  E  values  where  additional  values  of  Qtk(E),  not 
directly  found  in  the  table,  are  approximated  by  linear  inter¬ 
polation.  In  this  case  Tt^(z)  =  Tpg(z). 
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i.  for  t=36-55 : 


nnt(z)  =  nn3(z)  55 


(2t-71)  exp  [-  Yt/l^.Ttl(z)] 


55  r  Yt  -1 

Z  (2t-71)  exp  -  — 7~y 

t=36  L  KT*tlC  il 


with  Y  =  B2(t-36)  (t-35) 


j .  for  t=6-35 : 


Etl  =  B1  C(t-4)  (t'3)  '  (t'6)  (t'5)] 


k.  for  t=36-55 


Etl  =  B2  ^33) 


(t-36)  (t-35) 


where  0(X-X*)  = 


0  for  x  <  x 


1  for  x  >  x 


E-E. . „ 
r  ij 

alSO  =  J 
E1 


dE  D-.^E") 


,  q  A.  ,  E.  vPx;jA  /E  1/2  +  E.  .  \ 

"ith  v(E’E  j  ■  r/,2 


o*  T*Eijt 


[i 


_  fE  1/2  *  EU 


N. 

IJ* 


f°r  Eij£  <  E  <  E  and  E,  <  E1  <  E  -  E 
J  —  —  m  1  —  — 


=o  otherwise. 

The  integral  was  approximated  by  the  following  trapozoidal  type  logic. 
This  logic  was  also  used  wherever  the  speed  of  calculation  was  judged  to  be 
more  important  than  accuracy. 


Define  EMID^  = 


EMID„  = 


(ei+e2) 


2  2 


EMID,  = 


<E2+E3> 


3  2 


EMID, 


*-eneval-i  +  eneval) 


NEVAL 


emidneval+i  eneval 


Now  find  EMID..  ,  s.t 
N+l 


Then 


also 


and 


and 


and 


also 


EMIDm  <  E-E. . „  <  EMIDm  .. 
N  —  ljJl  N+l 


Qijl(E)  .  V(EMIDsll  -  EMIDs)  D..t(E,Es) 


*  CE-Eijt  -  EMIDn)  D.jt(E,EN) 


W“>E)  ■  >=ev(E)1/2  Qmj  (E) 


y  .  (z)  =  K  .  n.  (z) 
'eisv  J  ei  is 


f  (z,E)  =  K  n  (z) 
eev  '  ee  ev  ' 


ne(z)  =  Z  n.s(z) 
s 


G  .  .  (z ,E)  =  Z  G  .  .  (z ,E) 


PiJ 


pijm 
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with 


G  ..  (z,E)  =  n  .(z)  Q  ..  (E+E  ..  )*  I(z,E+E  ..  ) 
pijm  nj  xpijnr  pijm'  pijm' 

and  I(z,E)  =  I^CE)  exp  ^10^  secant  x  £  Qpaj  (E)  f  dz  nnj  Cz 
In  the  expression  for  Gp^(z,E),  (E) ,  QpijmCE)>  an(*  ICE)  are  bar 

fZm  11 

graphs.  Thus,  if  we  can  evaluate  J  dz  nnj (z  )  we  can  analytically 


r 


m 


evaluate  E.  E  G  ..(z,E)  since  G  ..(z,E)  is  piecewise  constant.  F . (z) 
i  j  piD  Pil  3 

f  zm  1  1 

/  dz1  n  .  (z1)  is  approximated  as  follows: 

» 7  a3 


The  logic  reads  in  a  set  of  z^  values  s.t.  z^<z^< . . ,<z^  . 

A  solution  is  found  at  one  of  these  z  values,  say  z..  Then: 

F . (z  )  =  0 
j  nr 

for  all  other  z . 1 s 

l 


F.(z.)  =  E  w_.  .  n  .(z.) 
3  1  t=i  t,i  nj v 


where  w.  .  are  the  weights  found  by  dividing  the  integral  into  sub-integrals 
usually  of  ORDINT  points.  Thus  an  ORDINT  point  integration  rule  is  used 
for  the  points  m,m-l , . . . ,m-0RDINT+l  and  then  points  m-ORDINT+1 ,m-0RDINT, . . . , 
m-20RDINT+2 ,  etc.,  until  m-k(ORDINT-l)  <  i  +  ORDINT  for  K  =  INTEGER  including 
0  at  this  point  the  remaining  sub-integral  is  evaluated  using  the  appropriate 
integration  rule  for  the  remaining  values. 


In  the  above  logic  n^(z)  and  Te(z)  are  guessed  values  and  all  values 
not  otherwise  described  are  input  parameters  or  are  described  by  a  table  of 
values  versus  E  or  z  whichever  is  applicable. 
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After  making  an  initial  guess  the  iteration  procedure  is  begun.  For 
each  iteration  the  values  of  ERR?^  and  the  C  matrix  are  calculated.  The 
procedure  first  calculates  ERR^  and  the  first  row  of  the  C  matrix.  This 
corresponds  to  the  values  due  to  the  approximation  to  the  equation 


SCALE 1  (  J  M  dE(E)1/2  H(z,E))/ne(z)  =  SCALE1 


where  SCALE1  is  an  input  parameter  adjusted  so  that  ERR^  is  the  same  order 
of  magnitude  as  the  other  ERR^'s.  The  integral 

•E 


/ 


dE(E)1/2  H(z,E) 


is  approximated  as 


M 


,1/2 


Z  (Ei)  '  (EMID.+1-EMID.)  H(z,E.) 
i=l 


and  therefore: 


M 


,1/2 


ERR,  =  SCALE1 *  (Z  (E.)  '  (EMID.  . -EMID. )  H(z,E.)/n  (z) 

1  .,1  1+1  1  1C 

1=1 


-  SCALE1 


dERR  M  n 

The  C..  values  =  *  *  =  SCALE1*(  Z  (E. )  '  (EMID  -EMID  )/n 

11  drl^Z  i_i  1  i+i  i 

-  (Ei)1/2  (EMID.+1-EMID.)  H(z,E.)/ne2  where  we  calculate 

3ne  ^nis(z) 

n  =  In.  (z)  and  tth-t — s-r-  =  £  — ,w_  -  as  follows: 

Q  -t  C  V  ' 


e  is 


3H(z,E.)  3H(z,E.) 


/ 


m 


ni2  =  %U  +  %ni2>  KGV  {  dE  EQr2^«  H(2>E)  +  [K211(Z) 


-1 


K24iU)]  nnl  +  tK233^Z^  +  K243*-Z^  "n3  +  K244  y  nn3 
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where  we  approximate 


r 


m 


dE  E  Qrs(E)  ll(z  ,E) 


as 


Q  M 
xrs 


2-n 


2^n~~  .£  (EMID.+1-EMID.)  (E.)  rS  H(z,E.) 


rs  i=l 


-n 


with  Qrs  and  nr<_  input  parameters  s.t.  Q^c(E)  =  Q^CE)  and 


rs 


rs 


'rs 


qpI2  "  qpI2(z)  +  f  dE  Gpil2(z,E) 


1 


with 


PiJ 


„(.)  -  / 


m 


dE  G  . . (z,E) 
PiJ 


J1 


r  m  rL 

where  /  dE  G  .  _(z,E)  and  /  m  dE  G  ,.(z,E)  are  approximated 

-/C  P-L1Z  j  pij 

1  bl 

in  the  explanation  of  the  initial  guess  and 

.E 

'E, 


as 


with 


r  e 

q  .  .  =  /  m  dE  G  . . (z,E) 
nemj  J  enij  *■  ’  J 


G  . . (z , E)  =  G  .  . (z,E)  -  G  .  . (z,E) 
enij  7  empj  J  emmj v  *  ’ 


where 


WZ-E)  *  Kev  VZ>  2  j l  J, dElElLijt(E 

ij£ 

G0„imJz.E)  =  Kev  n  .  (z)  EH(z,E)  Z  0(E-E..„)  Q..„ 
emmj v  nj v  '  a  13  V  *. 


Lijt<El'«  ’  2  'ijt  ta"'1 


'ljt'  '  1JJ 

1 


>*(E/Iijt)2  l.((EI-Eijt-E)/I..l]2 


described 


1  »E)  H(z jE1) 
(E) 


Further 


for  and  EjW-E. jt 

LijA(E^»E)  *  o  otherwise 

Here  Q^j^(E)  is  approximated  as  described  in  the  explanation  of  the 
initial  guess. 


E+E . 
ljJl 


is  approximated  as  follows: 


/  m  dE  e1  Liu<El,E)  HCz.E1) 


Let  EMID  <  E  +  E. . „  <  EMID  , 
t  —  ij£  t+1 


Then 


X 


dE1  E1  L..  (E1,E)  HCz.E1) 

E+E.. 
iJ  «- 


2 


M 


it.!  ‘EMIDs.rEMIly  Es  Eij£<Es’E>  hc*.es) 


and 


t  (EMIDMI-E-E1Jt)  Et  L  (E^.E)  HU.Ep 


/ 


dE  G  . . (z,E) 

F  eniJ 

1 


is  approximated  by 


M 


^  CENIDttl-EMIDt)  Genij(,.Et) 


This  leads  to  a  form 


M 

q  •  •  =  l  U  .  H(z,E  ) 
Menij  tj  v 


where  the  's  are  constants.  We  also  have 
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K  f  zl  =  CK 
211UJ  LK211  tT  (z)J 

n 


v  -  rv  ‘  1000^ 

^241  Z  LK241  (T  (z) 


K211(Z) 


v  c-i  -  rv  rJMLi 
K233LZj  233  lT  (z)J 
n 


if  r_'>  _  rK  rJOOO^ 
243UJ  "  LK243  lT  (z)J 
n 


1000 

K244^Z)  =  CK244  (T  (z)-* 
n 

Wlth  nn3(-Z')’  Y’  nnl(z')’  KeV’  Eijl’  *ijr  nrs’  CK211’  N211’  CK241’  CK233, 

n233 ’  C1C243’  n243 ’  CK244’  n244  inPUt  Pa™*ers- 
This  implies  n^2(z)  is  of  the  form 

M  /  M  _1 

ni2(z)  (C0N12  +  2  vt2  H(z,Et))  (  2  wt2  H(z,Et)  +  CON^) 

t=  1  \  t= 1 

8  ni2  /  M 

where  CON^,  vt2>  w^,  CON22  are  constants  and  =  v_2  /  2  *t2  H(z,Et) 

i  V  t  ^ 


+  CON, 


2)  ‘  • " 


i2‘CON12  .  ^  vt2  HU,Et))l  ^  »t2  H(Z,E 


t>  *  CON22  | 


nil<z)  *  (VlU)  *  %nil  *  K211(Z) 


*  K142(z)  nn2  *  K133(z)  n„3  *  1WZ)  V 


"nl  ni2>  /Kev  l  d 

'  IWZ)  *  "n 3V1 


dE  E  Qrl(E)  H(z ,E) 


if  -  rif 

K142(Z)  ~  CK144  T  (z)5 
n 


44 


W2) 


ck  r  iQQQ-v 
LK133  4’  (z y 


K142 ( Z) 


(_Z50_) 

4  (z)J 

nv  J 


_Z50_) 

4n(z)) 


for  Tn(z)  <  750 


for  (z)  750 


where  CK^,  n144,  CK^,  n^,  CK^,  n142a>  n142b  are  input  parameters. 


Further 


4il^  %il  ~  f 


f  dE  Gpil2(Z’E) 

4 


which  is  approximated  exactly  as  C z) . 
The  approximation  used  for 


r  E 

J  m  dE  E  Qrl(E)  H(z,E) 


has  been  previously  introduced.  With  these  approximations  n^  is  of  the 
form 


nil (z)  =  (C0Nn  +  I  vtl  H(z,Et)  +  CON21n.2(z)) 


k  i 


"tl  *  C0N31 


where  CON^,  v  ,  C0N21>  w  CON.^  are  constants  and 

9n  (z)  8n._(z)  / M  \ 

3H(z,E.)  (uil  +  C0N21  9H(z,F..)^  /  ^  Wtl  ll(z,Et) 
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Next 


>i  J  1  -  "u1 


M 


+  CON,,  l  -  w. , (CON, ,  +  Z  vtl  H(z,EJ  +  CON,,  n.,(z)) 

t=l 


M 

Z  w  H(z,E  )  +  CON 
t=l  * 


-2 


ni30)  -  (qpi3  +  qeni3  +  Ki33^z^  nn3  nn(z)  +  K233^ 


(Kev  / 
\  E1 


nn3  ni2(z)  (  Kev  J  m  dE  E  <WE)  H^Z'E^  +  K342(Z)  nn2 


with 


+  K344(Z)  y  nn3 


n 


K342^Z-*  “  CK344  (z)-1 

n 


342 


if  r-i  _  rv  '  1000' 
K344  Z  CK344  (T  (z) 


344 


where  n342’  E^344’  n344  are  ^nPut  parameters.  Using  approximations 

similar  to  those  used  for  n.j(z)  and  n^2(z)  this  equation  has  the  form. 

ni3(z)  -  (CON13  +  ^  vt3  H(z,Et)  ♦  CON23  n.^z) 


M 

+  CON33  ni2(z))f  tE1  Wt3  H(Z’Et} 


*  CQK4-^j 


-1 


where  CON^,  vt3,  CON23>  CON33>  wt3,  CON43  are  constants  and 
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3  n. 
i3 

3  H(z,E.) 


9nil(z)  (z) 

(vi3  +  C0N23  3H(z,E.)  +  CON33  aHCz.E.)3 


M 


-I 


Wt3  H(Z’V  +  C0N43  ]  -  Wi3(CON13 


M 

+  ^  Ut3  HC2*Et>  +  C0N23  "il^  +  C0N33  "i2<z” 


M 

Wt3  H(z’Et)  +  C0N43 


-2 


Next 


(■ 


ni4(z)  =  I  n.9  +  Ki aaW  y  n„J 


142''  ^  n2  144 1  "n3J  “il 

+  [K241(z)  nnl  +  K243(Z)  nn3  +  K244(z)  Y  ^  ni2 

E 

*  <K342(Z)  n„2  +  K344(z)  *  "nS1  ni3) /Kev  L  "  dE  E  <WE> 


H(z,E) 


Using  the  approximation 


S'  dE  E  <WE> 


this  equation  has  the  form 

n.4(z)  =  (CON14  n.^z)  ♦  CON24  n.2(z)  +  CON34  n.3(z)) 

M  _1 

Z  w  H(z,E  )) 
t=l  X  r 

where  CON^,  CON24,  CON34,  w  4  are  constants  and 


47 


3  n 


i4 


3  H(z,E.) 


3  n  .  (z)  3  n  (z) 

(CON,  „  „  ^  ■  +  CON 


14  3  H(z,E.) 


24  3  H(z ,E. ) 


3  n  (z)  /  M  \  -1 

*  C0N34  FTTUTeT))  »t4  H<z-M  wi4 

(CON14  nil(z)  +  CON24  ni2(z)  +  CON34  ni3(z)) 

M  '  "2 

E  w  H(z,E  ) 
t=l  z 


and  since 


3ne(z)  4  3n.s(z) 


ne(z)  E  n.^z)  5  3H(z,E  )  3H(z,E.) 

s=l  v  ’  t  s=l  1 


we  can  calculate 


3  ERR, 


3H(z,E.) 


which,  in  general  is  unequal  to  0  except  for  Ej=0. 


So  the  first  row  of  the  C  matrix  is  filled  except  for  the  first  element 
(at  this  point) 


C  =  /  0  X  X  X 


XXX 


/ 


,old 


After  the  first  row  of  the  C  matrix  and  ERR^  have  been  calculated, 

the  logic  calculates  the  second  row  of  the  C  matrix  and  ERR°lc*. 
be  followed  by  the  third  row  of  the  C  matrix  and  ERR°^,  etc.  The  row 


This  will 


of  the  C  matrix  and  ERR. 

J 


old 


uses  the  differential  equation  approximated  by 


a  difference  equation  for  E=E^  .  Thus  equations  for  E^.E^, . . .  are 

considered.  To  calculate  each  row  of  C  and  its  associated  ERR°^  value, 
the  effects  of  TERM1 ,  TERM2 ,  TERM3,  TERM4 ,  TERM5 ,  TERM6 ,  TERM11,  TERM12 , 
TERM13  and  TERM14  are  analyzed  separately.  The  effects  TERM7  with  TERM8 
and  TERM9  with  TERM10  on  each  row  of  C  and  its  associated  ERRol<*  value  are 
considered  in  pairs. 
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The  effect  of 


('2m  \ 

- ®_A 

m  .+m  / 
k  nj  e/ 


TERM1  =  E  I - -^-1  Ymn(z,E)L  3/2  H(z,E) 


on  the  C  matrix  and  ERR°*d  is  easily  calculated.  Here  v  (z,E}  =  n  fz) 
1/2  mn^  nj^ 

Kev  (E)  QjujjCE)  and  me>  ,  me>  n^  (z)  ,  and  Kev  are  input  parameters. 

Qmn CE)  is  a  table  of  values  versus  F.  which  is  linearly  interpolated 

for  any  missing  value.  The  effect  on  ERR?ld,  2  <_  i  £  NEVAL,  is  to  add 

the  factor 


J  \  n]  e 

and  the  effect  on  the  i^*1  row  of  C  is 


/2m  \ 

'  ?  )  y  (z,E.)  (E.)3/2  H(z,E.) 

mn  1  1  v  ,  ±J 


?  Wz-Ei>  <Ei>5/2  cit  s<‘-« 


Thus:  TERM1  adds  to  the  following  elements  of  the  C  matrix. 


0,0...,0,...,0 
0,X,0...,0,...,0 
0,0,X,0. . . ,0, . . . ,0 
0,0,0,X,0,. .. ,0, ... ,0 

o,o,...,o,...,o,x 

2  m 


The  effect  of  TERM2  =  E 

m.  +m  eis 
s  is  e 


e  1 

Y„; _ (z)  H(z,E)  on  the  C  matrix  and 


ERR0  is  also  calculated.  Thus  with  v* .  (z)  =  K  .  n.  (z)  and  m 

CIS  ° 1  i  c  v  '  /\ ' 


ei  is 


.  m.  ,  and 
e  is 


Kgi  as  input  parameters,  the  effect  of  TERM2  on  ERR^  with  2  <_  j  <  NEVAL,  is 


to  add  the  factor  E 
s 

j  row  of  C  is  found  as  follows: 


(1M  i 

Im.  +m  /  Yeis 
V  IS  e/ 


(z)  H(z, E) .  The  effect  of  TERM2  on  the 


*  3 (TERM2  for  E=E  ) 

add  to  C.t  -  C>t  3H(t,EE)  E_ 


(]i  "e  \  1 

(m.  +m  7  Yeis 
\  is  e/ 


(z)  6 (j  -i.)  +  l 


/2m  \ 

( - “  J 

Im.  +m  / 
\is  e/ 


*^mis^z^  ('ll  (2) 

K  .  v-  H  (z ,  E . )  =  TERM2t'  }  +  TERM21-  J  ; 

ei  3H(z,E£)  y 


here  the  effect  of  TERM2^  is  to  add  to  the  following  elements  of  the  C 


matrix 


0)0,0] . • . ,0, . . . ,0 

o,x,o, ... ,0, ...  ,0 

0,0, X,0, ... ,0, ... ,0 


0,0,0,. ..,o,...,o,x 


The  effects  of  TERM2^“'  is  to  add  the  C  matrix  wherever  any  of  the 
terms  3n^s/3H(z ,E^)  is  unequal  to  zero.  As  shown  in  the  explanation  of 
the  first  row  of  the  C  matrix  these  terms  are  generally  unequal  to  zero. 
Thus  the  effect  of  TERM2^  is  to  add  to  the  following  elements  of  the  C 


matrix. 


0,0,0, ... ,0, ... ,0 

x,x,x,...,x,...,x,x 

x,x,x,...,x,...,x,x 


x.x.x,. .. ,x,... ,x,x 


The  effect  of  TERM3  =  K  I  (z,E)  H(z,E)  on  the  C  matrix  and  ERR 

ee  o  ’  J 

is  also  calculated. 


rEn  i  i 1/2  i 

Thus  with  I  (z,E)  =J  dE  (E  )  H(z,E 


dE  (E  )  ll(z,E  )  which  is  approximated  as 


1=1 


S  (EMlDt+1-EMIDt)  (Et)1/2  H(z,Et)  +  (E^-EMID^)  (E£) 1/Z  H(z,E4) 

and  with  the  knowledge  that  Keg  is  an  input  parameter,  the  effect  of  TERM3 
on  ERR  with  2  <j  <  NEVAL,  is  to  add  the  factor  K  I  (z,E)  H(z,E.).  The 

J  .1  00  O  J 

effect  of  TERM3  on  the  j  row  of  C  is  found  as  follows: 

3 (TERM3  for  E=E.) 


1/2 


Add  to  C. „  =  C* „  = 

j*  n 


3H(z,E£) 


I  (z  ,E  .  )  6 (j -£) 
ee  o  j 


1/2 


Kee  (EMID5  +  1-EMID5)  (Ep )  '  H(z,E,)  for  £  <  j 


1+1 


,  1/2 


+  K  (E.-EMID.)  (E.)A/  H(z,E.) 

eevj  JJ  J 


for  £=j 
for  £  >  j 


=  TERN^1-*  +  TERM3*-2-1 


here  the  effect  of  TERM3^  is  to  add  to  the  following  elements  of  the  C 
matrix 

0,0,0, ... ,0, ... ,0 
0,X,0,.  . .  ,0,  .. .  .  ,0 
0,0, X,0, ... ,0, ... ,0 

o,o,o,x,o, . . . ,0, . . , 


0,0,0,. ..,...,o,...,o,x 

('21 

The  effects  of  TERM3V  ’  is  to  add  to  the  C  matrix  wherever  the  terms 
3  I  (z,E. ) 

'g  ^  are  unequal  to  zero  this  occurs  on  the  lower  diagonal 

except  for  the  first  row  and  the  first  column.  Thus  the  effects  of  TERM3^ 
is  to  add  to  the  following  elements  of  the  C  matrix 
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0,0,0, ... ,0, ... ,0 
0,X,0, ... ,0, ... ,0 
0,X,X,0, ... ,0, ... ,0 
0,X,X,X,0,...,0,...,0 


o,x,x,x,...,x,...,x,x 


The  effect  of  +TERM4 


„old 


.1 


dE  EG  ..(z,E)  is  to  add  "  constant  to 
E  P1J 


the  values  of  ERR^T  ,  for  2  <_  j  NEVAL.  This  constant  is  approximated  as 
described  in  the  method  for  determining  an  initial  guess  for  H(z,E). 


The  effects  of  -TERM5 


f  M  1  1 

-I  dE  E  G  ,.(z,E  )  is  described  in  finding 
J  enij v  ' 


,old 


n.  (z)  for  the  first  row  of  the  C  matrix.  The  addition  to  ERR.  ,  for 
1S  p  J 

/EM 

2  <  j  <  NEVAL,  is  equal  to  the  approximation  for 


~  J  dE1  I  G  . . (z jE1) . 


E. 

J 


3 


enij 


The  addition  to  the  C  matrix  does  not  include  the  first  row,  first  column 
and  last  row.  In  addition,  for  the  usual  input  data  some  of  the  first  and 
last  columns  have  no  contribution  from  this  term.  In  general,  however,  the 
addition  to  the  C  matrix  due  to  TERM5  has  the  following  form: 


0,0,0,. ...0,...,0 
o,x,x,x,. . . ,x, .. . ,x 
o,x,x,x,...,x,...,x 


o,x,x,x,. . . ,x, . . . ,x 

0,0,0,. .  .  ,0, _ ,0 


The  effects  of 


TERMO 


■/ 


M 


1/2 


dE1  Z  (E1)  Y  (z.E1)  H(z .E1) 
E  s 


■  ♦r 


i  i* 

(E  )  £n.s(z)  Kev  Qrs  (E)  rS 
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are  approximated  similarly  to  the  integral  described  in  the  first  row  of 
the  C  matrix.  So 


♦  J 


M  (E1)  E  n.s  (z)  Kev  Qrs(E)  "rs  HCz.E1) 


q  n.  M  2-n 

~  +  Kev  E  (-I® — — )  (  E  (EMID.  .-EMID.)  (E.  rS)  H(z,E.) 
«.  2-n  J  +  l  3  3  3 

b  rb  j=*+l 


+  <EMIIWV  (E^  °rS)  H(z, E^) 


which  is  added  to  ERR°^,  for  2  <_  l  <  NEVAL.  Analyzing  3ERR°^^/8H(z,Ej)  we 
see  that  TERM6  adds  terms  to  the  upper  triangular  portion  of  the  matrix  C, 
exluding  the  first  and  last  row. 

Thus  we  add  the  following  type  of  terms,  form  B,  to  the  C  matrix. 


0,0,0,. ..,0,...,0 
0,X,X,X,...,  ,X 
0,0,X,X,...,  ,X 


0,0,0, . . . ,0, . . . ,0,X,X 

0,0,0, . . . ,0, . . . ,0,0,0 


To  find  the  effects  of  TERM7  and  TERM8  we  note 
E. 


-TERM7 -TERM8 


t  k 


f  M  l  1  j 

=  -J  dE1  E  E  (EX  +  Etk)  y^Z.E^J 


tk' 


E  1/2 

H(z,E1  +  Etk)  ♦  J  M  dE1  E  E  (E1)  Y^U.E1)  Hfz.E1) 

fc  t  k 


fEM  11  i  i 

=  -  E  E  [  J  dE1  ftk(EA+Etk)  -  J  dE1  f^CE1)] 

t  k  E  E 


1  1  111  11 
where  f^CE1)  =  (E1)  y^U.E1)  H(z,EA)  =  E^n^z) *Kev*Qtk(EA) *H(z,EX) 


with  n  (z)  and  Kev  input  parameters  and  defined  in  the  description 

of  the  initial  guess.  We  approximate  -TERM7-TERM8  as  follows:  Use  an  input 
parameter  EREAD  and  rewrite  -TERM7-TERM8  as 


/-EREAD  ,  r  EREAD 

-0 (EREAD-E)  £  £  [J  dE1  -  J  dE1  f^E1)] 

t  k  E  ^E 


-EM/  dE1  ftkCE1*Etk) 

max (EREAD, E)  *^max(EREAD,E) 


f  dE‘  ftk<El 

•AiaxfF.RF.AD .  F.l 


) ]  where  we 


assume  EREAD-E  .  >  0  and  we  let  E,,-x=°. 

tK  M 

Now  approximate  f ( E1 )  by 


ftk(EHD  -  ftk<El> 

Etk 


and  cancel  terms  in  the  first  two  integrals  then  this  equation  can  be 
rewritten  as 


•EREAD  +  E 


tk 


-e(EREAD-E) 


t  k 


t  k 

“EREAD 

f  °°1 
tk1-  ' 

-  ftk 

dE  ftk^E 


11 -j; 


„E  +  E 


tk 


dE1  ftk(E1)] 


tk' 


We  wish  to  evaluate  this  expression  at  E^=E  ,  for  2  <_  i  <_  NEVAL,  and  with 

=AD  =  E  f 
n 

approximate 


EREAD  =  E  for  2  <  n  <  NEVAL  since  EREAD  is  so  set  by  the  logic  then  we  can 
n  — 


-6 (EREAD-E  )  £  £  [ 

t  k  E 


I 


E  +E. , 
n  tk 


dE1  f^tE1)]  by 


-8(E  -E()  E  E  U«n(E  .E  EH1D  )-En)  ftk(E  ) 
t  k 

NEVAL 

+  £  0 (EMID  -E  -E  ,  )  (min(E  +E^,  ,EMID  ,)-EMID  )f  .  (E  )] 

,  u  n  tk'  v  v  n  tk’  u+1  o  tkv  u' 1 

u=n+l 
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Since  f^k (Eu)  is  linearly  related  to  H(z,Eu)  then  derivatives  of  this 
approximation  with  H(z,Et),  for  1  <_  t  <_  NEVAL,  leads  to  additions  to  the 
C  matrix  for  row  2  through  row  n-1  on  the  ntE  column  and  possibly  for  a  few 
additional  terms  to  the  right  of  the  nc  column.  Similarly  we  approximate 


9  (E 


-EJ  *  E  / 

t  k  E„ 


VEtk 


dE1  f t k (E 1 )  by 


6<W  Z  l  [(1n(VEtk,EMID1+i)  -  Et)  ftk(Et) 

t  k 

NEVAL 

*  E  8,  tEMID  -E  -E  p  (min(E  .Etk,EMID  )  -  EMID  )  f  (E  )]  . 

U=£  +  l 

Since  f(E  )  is  linearly  related  to  H(z,Eu)  then  derivatives  of  this 
approximation  with  H(z,Et),  for  2  <_  t  <_  NEVAL,  leads  to  additions  to  the  C 
matrix  for  row  2  through  n-1  on  the  main  diagonal  and  possibly  a  few  columns 
to  the  right  of  the  main  diagonal.  Thus  we  add  to  ERR°lc^,  for  2  <_  j  NEVAL. 
The  value 

-0 (E  -E.)  E  Z  [(min(E  +E^.  ,EMID  J-E  )  f  ,  (E  ) 
n  y  .  1  v  n  tk’  n+1'  n ’  tkv  n' 

J  t  k 

NEVAL 

+  z  e(EMII)  -E  -E. , )  (min(E  +E  ,EMID  .)  -  EMID  )  f  ,  (E  )] 

,  v  u  n  tk'  v  v  n  tk’  u+1'  u'  tkv  u' J 

u=n+l 


+  6  (E  -E.)  Z  Z  [ (min(E . +E  ,,EMID .  ,)  -  E.)  f  ,  (E .) 
n  J  t  k  3  tk’  j+1  J  tk  3 


NEVAL 

+  z  e  (E.+E  k-EMID  )  (min(Eu+Etk,EMIDu+1)  -  EMIDu)  ftk CEu) ] 
u=j+l  J 

+  1  Z  Etk  ftkCmax(En’E j)) 
t  k  J 

where  f  ,  (E1)  =  E1  n  (z)  Kev  Q  .  (E1)  Hfz.E1)  and  we  add  to  the  C  matrix  a 
tk  nt  tk 

matrix  of  the  form 
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0,0,0,... 

0,X,  a  few  values, 0, 
0,0, X,  a  few  values. 


,0  X,  a  few  values, . 

,0  X,  a  few  values, . 


equation 


X,0,0, . . . 

o, x,o, . . . 
o,o,x,. . . 


0,0,0, . . . 


To  find  the  effects  of  TERM9  and  TERM10  we  note 


-TERM9-TERM10 


'  /  M  “e1  l  £  exp  ["  qT^U)]  <El*Etd 


Ytk(z,E1+Etk)  HCz.E1)  -  dE1  Z  Z 

E  t  k 


r  _Etk  1  fri,1/2 

oxp  m (E  1 


y  (z.H1)  6(E1-Etk)  HCz,E‘-Et|[) 


r  M  i  j  m 

=  z  X  [J  dE1  gtR(E1)  -  dE  gtk(E  -Etk) 


t  k  E 


where 


8(E  -Et|() 


8tk(E’)  '  W*' 

'(KT*Ttk'l))  ..1 


i  I  vTtk<2V  i 

,E1+E  , )  e  ^  7  HCz.E1) 


py 

vKt*t  k  (z)y  j  l  1 

V  tk  '  (E1+EtR)  nnt(z)  Kev*Qtk(E1+Etk)*H(z,E1) 


with  Kt,  Ttk(z),  nnt(z)  and  Kev  input  parameters  and  Qtk(E1)  defined  in  the 
description  of  the  initial  guess.  We  approximate  -TERM9 -TERM10  as  follows: 
Use  an  input  parameter  EREAD  and  rewrite  -TERM9-TERM10  as 

EREAD 

e(EREAD-E)  ZZ  [(/  dE1  g  (E1) 

tk  E  tK 


r 

-/ 


EREAD 


d£l  *tk  (El-Etk’  e‘El-Etk) 


/■*  OO  /**  OO 

V£  E  t  /  dE‘  stk<E‘>  -  dE‘  *tk  (E‘-Etk>l 

-max  (EREAD,  E)  J max (EREAD, E) 

where  we  assume  EREAD  -  E^k  >  0  and  we  let  E^-w.  Now  approximate  g^k(E^)  by 

Etk 

and  cancel  terms  in  the  first  two  integrals.  Now  this  equation  can  be 
rewritten  as 


« EREAD  -E 

0 (EREAD-E)  Z  Z  [  [  dE1  gtk(E1)  -  (  dE1  g(E1) 


t  k 


'EREAD-E 


tk 


E-E 


tk 


0 (E1) ]  +  Z  Z  Etk(gtk(“)  -  gtk (max (EREAD, E))) 
t  k 

but  gtk(“-)  =  ° 

We  wish  to  evaluate  this  expression  at  E=E  ,  for  2  <_  1  NEVAL,  and 
with  EREAD=En  for  2  n  <  NEVAL  (which  is  how  EREAD  is  set  by  the  logic), 
then  we  approximate 


0(E  -E  )  Z  Z  [ 
n  r  t  k 


ER-E 


dE1  gtk(E1)]  by 
tk 


+  9(EMIDu+l  •  En+Etk>  (EMIDu+l  ‘  "«ax(EMIDu,En-Etk))  g^)] 
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r 


Since  gt]c(Eu)  is  linearly  related  to  H(z,Eu)  then  derivatives  of  this 
approximation  with  Hfz.E^),  for  1  _<  t  £  NEVAL,  leads  to  additions  to  the  C 
matrix  for  row  2  through  row  n-1  on  the  nt^1  column  and  possibly  for  a  few 

|,L 

additional  terms  to  the  left  of  the  nL  1  column. 


Similarly,  we  approximate 


-e  (Er  e£)  £  £ 


dE1  gCE1)  0CE1)  by 
tk 


-9< W  9  l  l<Et-*“(EMIDu'Et-Etk» 

*  a(EMIDu.rVEtkK®IDu.l  -  "a^CBMID  Et-Etk))  g  (Ep] 

U=1 


Since  gtjc(Eu)  is  linearly  related  to  H(z,Eu)  then  derivatives  of  this 
approximation  with  H(z,Et),  for  1  _<  t  £  NEVAL,  leads  to  additions  to  the  C 
matrix  for  row  2  through  row  n-1  on  the  main  diagonal  and  possibly  a  few 
columns  to  the  left  of  the  main  diagonal. 

Thus  we  add  to  ERR?lc*,  for  2  <  j  <  NEVAL.  The  value 
3  ~  ~ 

0  (E  -E.)  £  £  [  (E  -max(E  -E. .  ,EMID  ))  g^,  (E  ) 
n  ,  1  n  ^ntk’  nJ J  6tkv  n' 

J  t  k 

n+1 

♦  E  e(EMIDu,rEn*Etk)  (EMIDutl  -  maxfEMID^.E^-E^p )  g^)] 


-6(En-E.)  E  E  [(E.-l»axCEMIDu,E.-Etk))  g^CEp 

*  V  6(EMID  -  E  -  E  )  (EHID  -  max(EMID  E  -E  ))  g  fc(E  )] 
11=  1  J  J  J 


+  f  l  Etk  gtktmax(En*Ej» 

t  k  J 


where 


J— tk  \ 

i  i  tVTtk<*V 

=  (E  +E^, )  e  N  1  tk  /  n 


*tk<E  )  =  <E  +Etk>  6 


KT*TtkCzV  1  1 

w  '  V(z)  Kev*Qtk (E  +Etk ) *H ( z * E  ) 
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and  we  add  to  the  C  matrix  a  matrix  of  the  form. 

0,0,0,  ...  >0,  ... 

X, X,0, . . . ,  a  few  values, X,0,  ... 

a  few  values, X,0, . . . ,a  few  values, X,0,  ... 
a  few  values, X, X, ... ,a  few  values, X,0,  ... 


0,0,0,. .. 
0,0,0,. . . 
0,0,0, . . . 
0,0,0, . . . 


,X,0,0,. . 
,0,X,0,.. 
,0,0,X,.. 
0,0,0, . . 


,0 

,0 

,0 

,0 


,0 

,0 

,0 

• • ,o,x 


th 

n  equation 


The  effects  of 


^  mp  T./2 

-TERM11  =  E  - - YmnO>E)  K^fz)  3H(z,E)/3E 


.  m  .+m  mni 
1  n3  e 


2  m 


=  E 


.  m  .+m 
1  nj  e 


„n.(z)  Kev(E)2  (E) 


are  approximated  by  approximating  at  E^=E  for  JQ  £  l  <  NEVAL.  Jq  is 

an  input  parameter  read  in  as  an  E  value  EjQ  l  <  EJ0  <_  E^q,  if  1  <  E^q 
<  NEVAL  and  if  EJ0=o  J0=1. 


/^H(z,E  )\  LL 

( -ir-4  =  y  Eu,* 

v  /aPP-  iM,l 


where  LL=min(°--D*--) ,  £-1  ,NEVAL-«.)  with  ORDINT  an  input  parameter  and  the 

E  set  to  be  the  constants  which  set  the  derivative  to  be  the  derivative 
u,  SL 

at  E^  of  the  lowest  order  polynomial  hitting  the  values  of  H(z,Efc+u)  for 
u  between  -LL  and  LL.  For  the  special  case 

■V 


Si  =  NEVAL 


/3H(Z,E  )\ 

( -  =  >  E  , NEVAL  H(z ,E  +u) 

’  \  3E  /  L—t,  u’  NEVAL 

'  app  u^l 
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If  «.<J  it  is  assumed  that  H(z,E)  can  be  approximated  in  the  region  around 

-( — - — ) 

VV 

by  Ke  with  K  a  constant  and  where 


M 


Te  =  (|)  (11605) 


f  dE  E3/2  H(z ,E)/n 


NEVAL 

(|)  (11605)  ^  E  EV2  (EMIDu+1-EMIDu)  H(z,Eu)  . 
e  u=  1 


In  this  case  we  approximate 


8H(z,E£) 

9E 


|3H(z,E£)^ 

\  3E  / 

app 


_C _ c _ .1  oo 

H(z,E  )  \*TeJ 

+e  E 


T  *KT 
e 


Ht2-W| 

\  E„+u  / 


u= -fcfc  -( 


VTe 


2  m 


Thus  we  add  to  ERR°ld  for  2  <  j  <  NEVAL  E  (-  --£-)  nnk(z)  Kev(E2)  Qmk 

J  If  n  If  o 


(E) 


/3H(z,E  )\ 
KTTn(‘Z')  l  9E  ) 

app 


k  nk  e 
If  0RDINT=7  and  EJO>0  the  terms 


/3H(z,E.)^ 
V  9  E 


_a£R 


3H(z,Es) 


lead  to  the  following  type  matrix  to  be  added  to  the  C  matrix. 
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J  = 
o 


0,0,0,. . . 

X,x,x, 

X,x,x, 


,0,... 

,X,.. . 

,x,... 


x,x,x,  ,x,.. 
0,0,0,. . .,o,x,x,x,x,x,x,x,o,. . 
0,0,0, . . . ,o,x,x,x,x,x,x,x,x,o, 
0,0,0,. .. ,o,o,x,x,x,x,x,x,x, . .. 


,0 

,X 

,x 


,X 

,0 

,0 

,0 


0,0,0,. . . 
0,0,0, . . . 
0,0,0, . . . 


0,X,X,X,X,X,X,X 

,o,x,x,x,x,x 

,o,x,x,x 


where  we  note 


(I  H(z,E4)^ 


aPP 


3H(z,Eo) 


for  i  <  J 


1  8Te  H(Z’E^ 

T7kT  6(EJl"Es)  +  3H(z ,E  ) 


e**T 


n 

+  E 
u=-JU 


-  ktt< 


3T 


2  Eu,*  3H(z ,E  ) 


+  e 


\  KT*Te  / 


E(z,E.  )  E  „  6(E  -E  )] 
1  t,+ir  u,£  r+u  s  1 


and  we  note 
3T 


3H(z,Es) 


(|)  (11605)  J-  Es3/2  (EMIDs+1-EMIDs) 
e 
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The  effects  of 


-TERM13  = 


2  Kee  3H(z,E) 

3  I2Cz»E)  — 3E - 


C  i  i  3' 1  1 

_ee  (  y  dE  (E  )  H(z,E) 


3H(z,E) 

3E 


r  i  i  3/2  i 

are  found  by  approximating  I  dE  (E  )  H(z,E  ) 


where  E  =  E  for  2  <  %  <  NEVAL  as  (I7(z,E  )  1  w  (E  )  '  H(z,E  )  where 

J6  A  A/  ,  U  .  J6  U  U 

app  u=l 

w  are  found  using  logic  similar  to  the  logic  used  in  finding  w  1  <  t  <  m 

U  j  Xi  l 


for  approximating 


F.  (z.)  f  m  dz1  n  .  (z) 
J  i  J  nj1- 


except  that  the  polynomials  start  at  E^  rather  than  z^ 
Using  these  approximations  we  add  to  ERR°^  for 


2  K  / 

2  <  j  <  NEVAL, — (I2(z,E  ))  (. 

•'  o  r~iTi  ' 


3H(z,E.)> 


-TERM13  at  E. 

The  value  of  — — •*-  for  2  <  1  <  NEVAL  and  1  <  U  <  NEVAL  are  found 
9Hlz,bul  3H(z,E.) 

by  chain  rule  differentiation.  For  0RDINT=7  the  derivatives  of  - v-g.--*1 — 


add  a  matrix  of  the  following  type  to  the  C  matrix. 


Finally,  the  effects  of 
2  K 


-TERM  14  =  — y^-  J(z,E)  3H(z  ,E)/3E  =  — y^-  E 


;3/2  /  dE1  H(z .E1) 


are  found  by  approximating 


e3/2  r  M  dEl 

J  C 


where  E=E^  for  2  <_  i.  <  NEVAL  as 


NEVAL 


(JCz.E^))  =  (E*)"'  (  E  Y  H(z.Eu)) 

app  u=  £ 


where  v  „  is  found  similarly  to  the  w.  used  in 
u,£  t 

a 

r  m  1  1  n 

F.-(z)  =  /  dz1  n  .(z1)  £  w.  n  .(z.) 

J  \  nj  t=i  t  njv 


-TERM14  at  E. 


The  value  of  — rrrT — ~  for  2  <  j  <  NEVAL  and  1  <  U  <  NEVAL 
9H(z,Eu)  -  -  - 


3H(z,E.) 

are  found  by  chain  rule  differentiation.  For  0RDINT=7  the  derivative  of  ( —  ■  ■  ) 


add  a  matrix  of  the  following  type  to  the  C  matrix. 


0,0,0, 

,0,0 

x,x,x, 

,x,x 

x,x,x. 

X 

X 

x,x,x. 

,x,x 

0,0,0, .. 

.  ,0,X,X,X,X,X,X,X,0  0,0 

0,0,0, 

,o,o,x,x,x,x,x,x,x,o,.,,,o,o 

0,0,0, 

,0,  ...0,X,X,X,X,X,X,X,0 

0,0,0, 

,0,  o,x,x,x,x,x,x,x 

0,0,0, 

,0,  o,x,x,x,x,x 

0,0,0, 

,0,  o,x,x,x 

0,0,0, 

,0,  0,0, 0,0 

^  1 
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